Identification of hub genes and potential ceRNA networks of diabetic cardiomyopathy

Diabetic cardiomyopathy (DCM), a common complication of diabetes, is defined as ventricular dysfunction in the absence of underlying heart disease. Noncoding RNAs (ncRNAs), including long noncoding RNAs (lncRNAs) and microRNAs (miRNAs), play a crucial role in the development of DCM. Weighted Gene Co-Expression Network Analysis (WGCNA) was used to identify key modules in DCM-related pathways. DCM-related miRNA-mRNA network and DCM-related ceRNA network were constructed by miRNA-seq to identify hub genes in these modules. We identified five hub genes that are associated with the onset of DCM, including Troponin C1 (Tnnc1), Phospholamban (Pln), Fatty acid binding proteins 3 (Fabp3), Popeye domain containing 2 (Popdc2), and Tripartite Motif-containing Protein 63 (Trim63). miRNAs that target the hub genes were mainly involved in TGF-β and Wnt signaling pathways. GO BP enrichment analysis found these miRNAs were involved in the signaling of TGF-β and glucose homeostasis. Q-PCR results found the gene expressions of Pln, Fabp3, Trim63, Tnnc1, and Popdc2 were significantly increased in DCM. Our study identified five hub genes (Tnnc1, Pln, Fabp3, Popdc2, Trim63) whose associated ceRNA networks are responsible for the onset of DCM.

www.nature.com/scientificreports/ is over-expressed in DCM, and silencing Kcnq1ot1 inhibits pyroptosis by regulating miR-214-3p and caspase-1 expressions. In addition, a recent study 15 showed that lncRNA, ZFAS1, acts as a ceRNA to sponge miR-150-5p and down-regulate CCND2, consequently promoting cardiomyocyte ferroptosis and DCM development. These studies indicate that noncoding RNA through the ceRNA hypothesis may play an important role in DCM.
In this study, we performed RNA-seq and miRNA sequencing on heart tissue from db/db mice to explore the transcriptome alterations in DCM pathogenesis, including lncRNA, miRNA, and mRNA. We further analyzed possible mechanisms by which altered immune cell infiltration may trigger DCM. This study provides a potential novel understanding of the pathogenesis of DCM and proposes several potential therapeutic targets. The research design is shown in Fig. 1. RNA-Seq pipeline. RNA was isolated from the left ventricle of three db/db mice and three db/m mice, respectively. RNA was used to generate barcoded complementary DNA (cDNA) libraries using the NEBNext Ultra RNA library prep with rRNA depletion (New England Biolabs). Indexed libraries were sequenced in 2 × 150 bp configuration on the Illumina HiSEQ platform. Briefly, first-strand and second-strand cDNA were synthesized from the input RNA, and single primer isothermal amplification (SPIA) of the resultant cDNA was performed, followed by mechanical shearing of double-stranded cDNA (ds-cDNA) (Covaris E220 ultrasonicator; Covaris Inc, MA). Subsequently, end repair and A-tailing of sheared cDNA and construction of unique barcoded libraries by addition of adapters and PCR amplification (8 cycles) were conducted. The Agencourt AmPure XP bead (Beckmann Coulter) purified libraries were quantified by quantitative PCR (q-PCR) and size distribution was checked using the Agilent TapeStation 2200 system. The libraries were subjected to paired-end 150 bp sequencing on HiSeq 2500 sequencing system (Illumina). For all samples, including public data, paired-end reads were mapped to the mm 10 genome assembly using RNA-seq aligner, STAR, with default parameters. Duplicate and multi-mapped reads were removed with samtools 19 . Following alignment, reads were counted per gene and fragments, per kilobase per million (FPKM), with featureCounts tool and GENCODE m28 and NONCODE v6 annotation 20,21 . Gene counts were then normalized and analyzed for identified differential expression mRNAs (DE-mRNAs) and lncRNAs (DE-lncRNAs) using the Bioconductor package edger, which is operated with a > 1.5-fold change and P-value < 0.05 cut-offs 22 . FPKM was used for Weighted Gene Co-Expression Network Analysis (WGCNA). miRNA-seq pipeline. miRNA was isolated from the left ventricle (3 db/db mice and 3 db/m mice). The libraries were generated using NEBNext®Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA). After library preparation and pooling of different samples, the samples were subjected to Illumina NextSeq500 sequencing. Raw data with fastq format containing N, low quality, 5' adapter contaminants, ploy A, and lowquality reads were removed from the reads. Raw data without a 3'adapter or the insert tag was also removed from the reads as well. The small RNA tags from rRNA, tRNA, small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA) were detected by mapping with the Rfam database and further mapped to the reference genome, mm10, by Bowtie software 23 . Known miRNA was identified by mapping to the miRBase database with R package, miRDeep2 24 . MiRDeep2 was also used to count the reads numbers mapped to each miRNA. R package edgeR was used for differential expression analysis. A cut-off with fold change ≥ 1.5 and P-value < 0.05 was used for identifying differential expression miRNAs (DE-miRNAs).

WGCNA.
We used the WGCNA package (version 1.60) in R to find and combine highly correlated mRNAs into mRNA modules, which permitted an examination for correlation between each module and diabetic complication-prone 25 . Before using the WGCNA, we used surrogate variable analysis (sva) package and objectoriented microarray and proteomics analysis (oompaBase) package to remove the batch effect between different samples 26 . The power of β = 5 (scale-free R 2 = 0.8) was selected as the soft threshold to ensure a scale-free network ( Figure S1). The dynamic tree cutting method was used to cluster the mRNAs in layers, using 50 as a minimum size cutoff, and the cut height = 0.3 was applied to merge highly similar modules. Different mRNA modules were labeled with varying colors; the gray module represents mRNAs that could not be merged 27 . Pearson correlation analysis was implemented to evaluate the correlation between mRNAs and susceptibility to diabetic complications.
Functional enrichment analysis. Gene Ontology (GO), including cellular components, molecular function, biological process, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway 28,29 , were analyzed using the R package clusterProfiler (version 3.2.14), as described previously 30 . Gene set enrichment analysis (GSEA) was performed using KEGG pathway annotation data, analyzed with the package clusterProfiler, and displayed on ridgeline plot 31 . Cell senescence related genes were download from CSGene database 32  www.nature.com/scientificreports/ regeneration-related genes were download from Regeneration Roadmap 33 . DE-lncRNAs have a high positive correlation (correlation coefficient > 0.9) with cell senescence-related genes and heart regeneration-related genes, identifying potential cell senescence-related or heart regeneration-related lncRNAs.
Immune cell infiltration. The immune cell infiltration status was acquired from bulk RNA-sequencing data by applying the single-sample gene set enrichment approach to the transcriptomes based on singlesample GSEA (ssGSEA), and ultimately acquired through the R package, gene set variation analysis (GSVA) 34 (Table S1.). A total of 25 immune cells were used in this study, which include activated B cell, activated CD4 + T cell, activated CD8 + T cell, central memory CD4 + T cell, central memory CD8 + T cell, effector memory CD4 + T cell, effector memory CD8 + T cell, immature B cell, memory B cell, T follicular helper cell, Tγδ, Th1, Th2, Th17, Treg, natural killer (NK) cell, eosinophil, activated dendritic cell (DC), immature DC, neutrophil, plasmacytoid DC, macrophage, mast cell, monocyte, and NKT cell 35 .
Prediction of target miRNAs and construction of ceRNA networks. According to the ceRNA hypothesis, miRNAs can induce gene silencing and down-regulate gene expression by binding to mRNA, and lncRNAs that are rich in miRNA binding sites can act as a miRNA sponge, leading to changes in expression levels of miRNA-target mRNAs. We first identified differentially expressed miRNAs (DE-miRNAs) in DCM and then calculated the Pearson correlation coefficient between these DE-miRNAs and hub genes. DE-miRNAs paired with correlation coefficient < -0.5 and P < 0.05 were considered as potential miRNAs that can downregulate the expression of the target hub gene. The ability to bind target hub gene was established by RNAhybrid (energy threshold = -20, G: U in seed) 36,37 . Furthermore, DE-lncRNAs that have a high positive correlation coefficient (correlation coefficient > 0.99) with hub genes that bound the same miRNAs were used to construct ceRNA networks. The result was presented as a mRNA-miRNA co-expressed network and ceRNA networks by Cytoscape 38 .
qRT-PCR assay. Transcript One-Step gDNA Removal (YEASEN, shanghai, China) and cDNA Synthesis SuperMix (YEASEN, shanghai, China) were employed during the reverse transcription process according to the manufacturer's instructions. cDNA was amplified with SYBR q-PCR Master Mix (EnzyArtisan, shanghai, China) in StepOnePlus (Applied Biosystems) equipment. The 2-ΔΔCt method was used to calculate the relative lncRNA expression, and each hole was repeated three times to ensure quantitative accuracy. The sequences of primers used for qRT-PCR are listed in Table S2, and GAPDH was used as an internal reference gene.

Statistics analysis.
Statistical analyses were accomplished using R. The student's t-test was used for comparison between two independent groups. A P < 0.05 was considered statistically significant.
Ethics approval and consent to participate. All animal experiments were performed according to procedures approved by the Laboratory Animal Ethics Committee of Southwest Jiaotong University.

Gene expression analysis of DCM.
To distinguish expressed genes between DCM and healthy controls, we performed RNA-seq experiments using db/db mice model and its corresponding control db/m model. Our results reflect that 1493 genes were up-regulated and 1609 genes were down-regulated in DCM mice compared with control mice (Fold Change ≥ 1.5, P < 0.05) ( Fig. 2A and Table S3). We performed the KEGG pathway analysis to determine possible regulation mechanisms of these differentially expressed genes (DEGs) in the process of DCM. The top 10 pathways that involved up-regulated DEGs included oxidative phosphorylation and reactive oxygen species (Fig. 2B). On the other hand, the top 10 pathways associated with down-regulated DEGs included NOD-like receptor signaling pathway (Fig. 2C). Similar results were obtained by performing GSEA of KEGG pathway. The results reveal that the genes associated with oxidative phosphorylation and reactive oxygen species were mainly up-regulated in DCM mice (Fig. 2D). GO analysis found up-regulated DEGs were associated with some biological processes, which include ATP metabolic process, mitochondrion organization, and mitochondrial ATP synthesis, while down-regulated DEGs were associated with biological processes including the regulation of immune effector (Fig. 2E). Besides, the analysis also found abnormal expression in an abundant of genes associated with senescence and heart regeneration in db/db mice ( Figure S3). Since DCM can be facilitated by alterations in both adaptive and innate immune systems, we further examined the difference of immune cell infiltration in left ventricle between DCM mice model and control model. The results reveal that activated CD8 + T cells infiltration was significantly increased in the DCM compared to control, while the infiltration of memory B cells, NKT natural killer cells, monocytes, and mast cells was significantly decreased, indicating an impaired profile of innate immunity in DCM (Fig. 2F). www.nature.com/scientificreports/  The results demonstrate that up-regulated DEGs in DCM were enriched in multiple signaling pathways, including oxidative phosphorylation and reactive oxygen metabolism, as compared to other diabetic complications (Fig. 3A). In addition, we found that down-regulated DEGs in DCM were associated with NOD-like receptor signaling and chemokine signaling, while down-regulated DEGs in other diabetic complications mainly focused on other pathways, such as fluid shear stress, atherosclerosis, nucleocytoplasmic transport, and nucleotide excision repairment (Fig. 3B). Given the discovery that down-regulated DEGs in DCM is highly associated with the host innate immunity pathway, we further performed immune cell infiltration analysis to determine specific immune cell subtypes. The results show that the infiltration of activated CD8 + T cells was significantly increased in DCM, as compared to other diabetic complications (Fig. 3C). This indicates that excessive specific immune cell infiltration within the heart may contribute to the development of DCM.

The difference of KEGG pathway and immune cell infiltration between
Identification of key module promoting DCM development. Next, we performed WGCNA to further determine the key factors regulating DCM development. The sample cluster tree was constructed and demonstrated in Fig. 4A. Through WCNA analysis, 15 co-expression modules were constructed (Fig. 4B).
Module-trait relationship analysis revealed that purple modules were positively correlated to DCM, while the magenta module was negatively associated with DCM (Fig. 4C). In addition, gene expression in the purple module was significantly up-regulated in DCM compared to other diabetic complications. Contrarily, gene expression in the magenta module was significantly down-regulated in DCM compared to other diabetic complications (Fig. 4D). Altogether, these results suggest that genes in the purple and magenta module may affect DCM development.

Functional enrichment analysis of Co-DEGs.
To further explore the biological function of co-expressed DEGs (Co-DEGs), Co-DEGs were first obtained from the intersection of purple and magenta modules (Fig. 5A). Next, KEGG pathway analysis revealed that Co-DEGs in the purple module were mainly enriched in the cardiac muscle contraction pathway, while Co-DEGs in the magenta module were mainly enriched in extracellular matrix (ECM) receptor interaction pathway (Fig. 5B,C). GO Biological Process (BP) analysis found that the Co-DEGs in the purple module were mainly enriched in biological processes involved in muscle system process, and Co-DEGs in the magenta module were mainly enriched in biological processes in amino acid transportation (Fig. 5D,E).

Hub gene identification and validation.
In order to screen out the core genes in DCM, Co-DEGs in the intersection of purple and magenta modules were used to construct a co-expression network by WGCNA. The resulting data file was then processed with Cytoscape and demonstrated in Fig. 6A. In this largest connected master network, we identified 5 hub genes by cytoHubba (Table S4). Afterward, we performed real-time PCR to check the expression of the five hub genes in the cardiomyocytes derived from either db/db or db/m mice. The results (Fig. 6B) show that there were significant increases in Phospholamban (Pln), Fatty acid binding proteins 3 (Fabp3), Tripartite Motif-containing Protein 63 (Trim63), Popeye domain containing 2 (Popdc2), and Troponin C1 (Tnnc1) gene expression in DCM compared to control (P ≤ 0.05). All these results were analyzed according to bioinformatics methods.

Prediction of target miRNAs based on hub genes. To further delineate the network between
identified hub genes and potential target miRNAs, we used miRNA-seq to predict the target miRNAs of the hub genes mentioned above. First, we identified 52 DE-miRNAs (Fold Change ≥ 1.5, P < 0.05) in DCM (Table S5). According to these results, we predicted that the expression of hub genes was most likely regulated by miRNAs. Under this prediction, we constructed a co-expressed network between 55 hub genes and miRNAs through Cytoscape (Fig. 7A). These hub genes were linked together by shared miRNAs. Next, we performed KEGG enrichment analysis to determine the functional role of these miRNAs in the pathogenesis of DCM. The results show that the target genes of these miRNAs were mainly involved in the TGF-β signaling pathway and Wnt signaling pathway (Fig. 7B). Moreover, GO BP enrichment analysis found that these miRNAs participated in the signaling of TGF-β and glucose homeostasis (Fig. 7C). Collectively, these results suggest that these miRNAs may promote the DCM process by affecting metabolic pathways and hub gene expressions.
Construction of CeRNA networks. Subsequently, according to the ceRNA hypothesis, we constructed ceRNA networks based on mRNA-miRNA co-expression network and related lncRNAs (Fig. 8). This network comprises 48 nodes and 517 edges and consists of five hub genes, 24 lncRNAs, and 15 miRNAs. Annotations of lncRNAs in ceRNA-network are shown in Table S6.  41,42 . We found abnormal expression in an abundant of cell senescence-related genes in DCM, indicating the close association and impact of senescence on DCM 43 .  www.nature.com/scientificreports/ Additionally, among 5 genes related to myocardial repair in mice (Hdac4, Mesp1, Ngf, Pim1, and Yap1), Ngf and Pim1 were significantly down-regulated, implying a decreased repair ability of cardiomyocytes in DCM and suggesting a potential treatment target to restore repairment in DCM.
In addition, we further performed immune cell infiltration analysis and revealed that activation of CD8 + T cells plays a predominant role in the development of DCM. Therefore, the development of diabetes is associated with abnormalities in the immune system. However, the role of immune cells in diabetic myocarditis is still unclear. T lymphocyte infiltration into the myocardium has been observed in a left coronary artery occlusioninduced murine myocardial infarction model and a transverse aortic constriction-induced murine pressure overload model 44 . Tang et al. 45 reported that CD8 + T cells in ischemic failing human hearts may contribute  www.nature.com/scientificreports/ to the progression of heart failure. Abdullah et al. 46 showed that conditional T-cell sS1p1 knockout mice that exhibited sustained deficiency of both CD4 + and CD8 + T cells, had improved cardiac function and alleviated cardiac fibrosis after 11 weeks of diabetic induction, indicating that T cell Ss1p1 activation exacerbates fibrosis under hyperglycemia. Although current knowledge supports that CD4 + T cells play a more important role in the development of DCM, mainly via the subtype of CD4 + Foxp3 + T cells 47 , we speculate that CD8 + T cells may be involved in the development of DCM via the following mechanisms: (1) CD8 + T cell can directly damage cardiomyocytes via its cytotoxicity effect; (2) CD8 + T cell can regulate macrophage migration via stimulating the production of nitric oxide; (3) CD8 + T cell can up-regulate CD11b, CD64, and CD62L on neutrophils mainly through the secretion of inflammatory factors and resultantly maintain their survival. The current study focuses on individual diabetic complications. We analyzed the differences between tissues affected by diabetic complications by comparing them with public databases. As expected, there were differences in pathways and immune cell infiltration for each complication abnormality. We further clarified the key genes and key modules responsible for DCM by WGCNA, which were mostly differentially expressed in DCM and normally expressed in other diabetic complications. Enrichment analysis showed that up-regulation of genes related to DCM in key modules were mainly involved in calcium signaling, senescence pathways and hypertrophy related pathways, consistent with our previous comparison on DCM with DPN and DN. A recent study implementing RNA sequencing on STZ mice type 2 diabetes DCM model demonstrated that these pathways were down-regulated, which comments on the different underlying pathogenic mechanisms in type 2 diabetes, such as insulin resistance for db/db mice and islet β-cell reduction in STZ mice, that may influence DCM 48 .
Furthermore, we identified and validated 5 hub genes (Tnnc1, Pln, Fabp3, Popdc2, and Trim63) from the DCM-related key module. Mutations in Tnnc1, a complex that is known as Cardiac Troponin C and contains the component, troponin C, was confirmed to be associated with hypertrophic or dilated cardiomyopathy 49,50 . Pln is a 52-amino acid sarcoplasmic reticulum (SR) membrane protein expressed abundantly in cardiac muscle and a crucial regulator of cardiac function by modulating the rate of cardiac relaxation and size of the SR Ca2+ store 51,52 . Increased expression of Pln has been identified to reduce cardiac contractility and correlates with the overexpression of NF-SLN, raising the possibility that induced expression of SLN in human hearts can impair cardiac function 53,54 . Jia et al. 55 also reported that the expression of Pln was significantly increased in a time-dependent manner in diabetic groups. Fatty acid-binding protein 3 (Fabp3) participates in cell metabolism by binding free long-chain fatty acids (LCFAs) and transporting them for cell metabolism 56 . Fabp3-defect exacerbates cardiac hypertrophy and heart dysfunction, but over-expression of Fabp3 can up-regulate the phosphorylation of the MAPK signaling pathway and decrease phosphorylated Akt levels, which may account for the augmentation of apoptosis and remodeling after myocardial infarction 57,58 . Popdc2, one of the Popeye domain-containing (Popdc) gene families, was highly expressed particularly in the sinoatrial node of the mouse and represented as a novel arrhythmia gene for cardiac conduction disorders 59,60 . A recent study 61 showed that Popdc2 was a fasting-induced gene, which suggests that the abnormal expression of popdc2 may be related to blood glucose. Trim63, also known as MuRF1, was significantly increased not only in cardiac muscle of diabetic mice, but also in diabetic limb muscle and STZ-Diabetes 62,63 . Previous studies have demonstrated that abnormalities in these genes were strongly associated with the development of diabetes or cardiovascular disease. Our results classify these genes as key factors in the pathogenesis of DCM and potential drug targets for DCM treatment. However, there are only a few reports on the regulation of these genes with miRNA or ceRNA, except Trim63 64 .

Conclusion
Overall, we systematically analyzed the characteristics of both mRNA and noncoding expression profiles in DCM. We identified the potential mechanisms and the hub genes in DCM pathogenesis, namely Phospholamban (Pln), Fatty acid binding proteins 3 (Fabp3), Tripartite Motif-containing Protein 63 (Trim63), Popeye domain containing 2 (Popdc2), and Troponin C1 (Tnnc1). Additionally, we also found potential mRNA-miRNA and ceRNA networks in DCM. Our results shed light on further studies of DCM pathogenesis and on the discovery of DCM therapeutic targets.

Data availability
The datasets generated and/or analyzed during the current study are available in the GEO  www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.